function [cr,pr] = blsrho(so,x,r,t,sig,q) 
%BLSRHO Black-Scholes sensitivity to interest rate change. 
%   [CR,PR]= BLSRHO(SO,X,R,T,SIG,Q) returns the rate of change in value of
%   securities with respect to interest rates.  SO is the current security
%   price, X is the exercise or strike price, R is the interest rate, T is
%   the time until maturity expressed in years, SIG is the volatility (standard
%   deviation), and Q is the dividend rate.  The default Q is 0.  CR is the
%   call option rho and PR is the put option rho.
%       
%   Note: This function uses normcdf, the normal cumulative
%         distribution function in the Statistics Toolbox. 
% 
%   For example, [c,p] = blsrho(50,50,.12,.25,.3,0) returns 
%   c = 6.6686 and p = -5.4619. 
% 
%   See also BLSPRICE, BLSDELTA, BLSGAMMA, BLSTHETA, BLSVEGA, BLSLAMBDA. 
% 
%   Reference: Hull, Options, Futures, and Other Derivative Securities, 2nd
%              edition, Chapter 13. 
 
%       Copyright 1995-2006 The MathWorks, Inc.
%       $Revision: 1.6.2.4 $   $Date: 2009/04/15 23:07:21 $ 
 
 
if nargin < 5 
  error('finance:blsrho:missingInputs',sprintf('Missing one of SO, X, R, T, and SIG.')) 
end 
[m,n] = size(so); 
if nargin < 6 
  q = zeros(m,n);   % default dividend rate 
end 
if any(so <= 0 | x <= 0 | r < 0 | t <=0 | sig < 0) 
  error('finance:blsrho:invalidInputs',sprintf('Enter SO, X, and T > 0. Enter R and S >= 0.')) 
end 
 
message = blscheck('blsrho', so, x, r, t, sig, q);
error(message);


% Perform scalar expansion & guarantee conforming arrays.
try
    [so, x, r, t, sig, q] = finargsz('scalar', so, x, r, t, sig, q);
catch
    error('Finance:blsrho:InconsistentDimensions', ...
        'Inputs must be scalars or conforming matrices.')
end

% blspriceeng works with columns. Get sizes, turn to columns, run engine,
% and finally turn to arrays again:
[m, n] = size(so);

% Double up on fcn calls since blsprice calculates both calls and puts. Do
% this only if nargout>1
NumOpt = numel(so);
callSpec = {'call'};
callSpec = callSpec(ones(NumOpt,1));
putSpec = {};
if(nargout)>1
    putSpec = {'put'};
    putSpec = putSpec(ones(NumOpt,1));

    % double up the rest of the input args
    [so, x, r, t, sig, q] = deal([so(:);so(:)], [x(:);x(:)], [r(:);r(:)], ...
        [t(:);t(:)], [sig(:);sig(:)], [q(:);q(:)]);
end

OptSpec = [callSpec;putSpec];
OutSpec = {'rho'};

% call eng fuction
rho = blspriceeng(OutSpec, OptSpec, so, x, r, t, sig, q);

% Now separate calls from puts
cr=reshape(rho{1}(1:NumOpt), m, n);
if(nargout>1)
    pr = reshape(rho{1}(NumOpt+1:end), m, n);
end